All data aggregated but two different fundamental assumption sets depending on how constraints associated with unused points. Unused can be fully randomly sampled, i.e. could in theory overlap with potential animal relocation at a different point in space or time. Alternatively, unused points sampled from environment can be constrainted to ensure that they never overlap with another relocation instance.
#data with replace means can overlap (same assumption as work in Carrizo)
#data <- read_csv("data/data_with_replace.csv")
data <- read_csv("data/tidy_data.csv")
data <- data %>%
select(rep, status, year, ground, mesohabitat, Aspect, Slope, Elev, KDF, Shrub_Cov, NDVI, Solar, aspect.cat)
data$year <- as.character(data$year)
data
## # A tibble: 6,708 x 13
## rep status year ground mesohabitat Aspect Slope Elev KDF
## <dbl> <dbl> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 2 0 2016 above open 217. 3.58 708 3.24e-3
## 2 3 0 2018 below open 218. 4.08 709 3.77e-3
## 3 7 0 2018 above shrub 247. 2.73 697 7.69e-5
## 4 8 0 2018 above shrub 202. 3.85 710 3.72e-3
## 5 9 0 2018 above shrub 225 0.506 698 3.26e-5
## 6 11 0 2016 below open 225 3.54 709 3.62e-3
## 7 16 0 2018 above shrub 288. 1.13 692 6.90e-6
## 8 21 0 2016 above open 217. 3.58 708 3.24e-3
## 9 24 0 2017 above shrub 63.4 1.60 698 4.10e-6
## 10 27 0 2018 above shrub 297. 3.20 698 0.
## # ... with 6,698 more rows, and 4 more variables: Shrub_Cov <dbl>,
## # NDVI <dbl>, Solar <dbl>, aspect.cat <chr>
#function
#global availability (m=0) and bootstrap (B=99)
# m = 0 indicates sample from anywhere in dataframe versus structure or blocks
model.rspf <- function(x) {rspf(status ~ Shrub_Cov + NDVI + Slope + Elev + Aspect + Solar, x, m=0, B = 99, model = TRUE)
}
m7 <- model.rspf(data)
summary(m7)
##
## Call:
## rspf(formula = status ~ Shrub_Cov + NDVI + Slope + Elev + Aspect +
## Solar, data = x, m = 0, B = 99, model = TRUE)
##
## Resource Selection Probability Function (Logistic RSPF) model
## Non-matched Used-Available design
## Maximum Likelihood estimates
## with Nonparametric Bootstrap standard errors (B = 99)
##
## Fitted probabilities:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.8263 0.9963 0.8236 0.9997 1.0000
##
## Coefficients (logit link):
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.228e+01 5.510e-01 40.440 < 2e-16 ***
## Shrub_Cov 9.580e+00 3.521e-01 27.205 < 2e-16 ***
## NDVI 9.132e+00 7.249e-01 12.597 < 2e-16 ***
## Slope 2.428e-01 3.575e-02 6.792 1.1e-11 ***
## Elev -2.663e-01 1.322e-02 -20.138 < 2e-16 ***
## Aspect -3.043e-04 1.666e-03 -0.183 0.855
## Solar 4.469e-04 2.564e-05 17.429 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Log-likelihood: -2.647e+04
## BIC = 5.3e+04
##
## Hosmer and Lemeshow goodness of fit (GOF) test:
## X-squared = 292.4, df = 8, p-value < 2.2e-16
plot(m7)
mep(m7) # marginal effects similar to plot but with CIs
kdepairs(m7) # 2D kernel density estimates
#split-map to explore key levels####
#year####
m.year <- data %>%
split(.$year) %>%
map(~ model.rspf(.))
m.year
## $`2016`
##
## Call:
## rspf(formula = status ~ Shrub_Cov + NDVI + Slope + Elev + Aspect +
## Solar, data = x, m = 0, B = 99, model = TRUE)
##
## Resource Selection Probability Function (Logistic RSPF) model
## Non-matched Used-Available design
## Maximum Likelihood estimates
##
## Coefficients (logit link):
## (Intercept) Shrub_Cov NDVI Slope Elev
## -2.059e+00 1.653e+01 4.346e+00 -6.792e-02 1.188e-02
## Aspect Solar
## -1.447e-03 -2.424e-05
##
##
## $`2017`
##
## Call:
## rspf(formula = status ~ Shrub_Cov + NDVI + Slope + Elev + Aspect +
## Solar, data = x, m = 0, B = 99, model = TRUE)
##
## Resource Selection Probability Function (Logistic RSPF) model
## Non-matched Used-Available design
## Maximum Likelihood estimates
##
## Coefficients (logit link):
## (Intercept) Shrub_Cov NDVI Slope Elev
## -5.921e+01 1.347e+00 2.099e+01 -8.013e-02 -6.508e-03
## Aspect Solar
## 1.274e-03 1.497e-04
##
##
## $`2018`
##
## Call:
## rspf(formula = status ~ Shrub_Cov + NDVI + Slope + Elev + Aspect +
## Solar, data = x, m = 0, B = 99, model = TRUE)
##
## Resource Selection Probability Function (Logistic RSPF) model
## Non-matched Used-Available design
## Maximum Likelihood estimates
##
## Coefficients (logit link):
## (Intercept) Shrub_Cov NDVI Slope Elev
## 6.967e+01 2.933e+00 1.399e+01 7.802e-02 -2.341e-01
## Aspect Solar
## 3.633e-05 2.521e-04
#ground####
m.ground <- data %>%
split(.$ground) %>%
map(~ model.rspf(.))
m.ground
## $above
##
## Call:
## rspf(formula = status ~ Shrub_Cov + NDVI + Slope + Elev + Aspect +
## Solar, data = x, m = 0, B = 99, model = TRUE)
##
## Resource Selection Probability Function (Logistic RSPF) model
## Non-matched Used-Available design
## Maximum Likelihood estimates
##
## Coefficients (logit link):
## (Intercept) Shrub_Cov NDVI Slope Elev
## 10.2942123 10.2087814 5.4884795 0.3357263 -0.3496654
## Aspect Solar
## 0.0005631 0.0006408
##
##
## $below
##
## Call:
## rspf(formula = status ~ Shrub_Cov + NDVI + Slope + Elev + Aspect +
## Solar, data = x, m = 0, B = 99, model = TRUE)
##
## Resource Selection Probability Function (Logistic RSPF) model
## Non-matched Used-Available design
## Maximum Likelihood estimates
##
## Coefficients (logit link):
## (Intercept) Shrub_Cov NDVI Slope Elev
## 24.6215758 8.3778257 10.8635911 0.1259252 -0.2341786
## Aspect Solar
## -0.0003167 0.0003786
#mesohabitat
m.meso <- data %>%
split(.$mesohabitat) %>%
map(~ model.rspf(.))
m.meso
## $open
##
## Call:
## rspf(formula = status ~ Shrub_Cov + NDVI + Slope + Elev + Aspect +
## Solar, data = x, m = 0, B = 99, model = TRUE)
##
## Resource Selection Probability Function (Logistic RSPF) model
## Non-matched Used-Available design
## Maximum Likelihood estimates
##
## Coefficients (logit link):
## (Intercept) Shrub_Cov NDVI Slope Elev
## 36.6829900 -0.4861748 3.6675175 0.2149919 -0.2715408
## Aspect Solar
## 0.0018867 0.0004218
##
##
## $shrub
##
## Call:
## rspf(formula = status ~ Shrub_Cov + NDVI + Slope + Elev + Aspect +
## Solar, data = x, m = 0, B = 99, model = TRUE)
##
## Resource Selection Probability Function (Logistic RSPF) model
## Non-matched Used-Available design
## Maximum Likelihood estimates
##
## Coefficients (logit link):
## (Intercept) Shrub_Cov NDVI Slope Elev
## -2.843e+00 2.377e+01 8.713e+00 -8.904e-02 -1.748e-02
## Aspect Solar
## 2.074e-03 2.966e-05
A few instances in literature calculate log odds ratio to estimate effect sizes. Contrast of used by unused probabilities.
#build dataframe
#model <- as_tibble(m7$model)
#fitted_values <- as_tibble(m7$fitted.values)
#odds <- bind_cols(model, fitted_values)
#odds
#write_csv(odds, "data/odds.csv")
#library(esc)
#effect_sizes(odds, t = value, fun = "esc_bin_prop", es.type = "or")
#just to run this section only
library(tidyverse)
#library(scales)
#odds <- read_csv("data/odds.csv")
#log odds ratios####
#log(p(A),p(not-A)) typically for A as event or cell frequency
#log((p1/(1-p1))/(p2/(1-p2))) for proportion successful in each group
odds_df <- read_csv("data/log_odds.csv") %>%
select(-aspect)
odds_df_long <- odds_df %>%
gather(factor, measure, 1:5)
ggplot(odds_df_long, aes(measure, log_odds)) +
geom_smooth(method = lm, color = "darkgreen") +
facet_wrap(~factor, scales = "free") +
labs(x = "", y = "log odds ratio")
library(broom)
m <- odds_df_long %>%
group_by(factor) %>%
do(fit = lm(log_odds~measure, data =.))
tidy(m, fit)
## # A tibble: 10 x 6
## # Groups: factor [5]
## factor term estimate std.error statistic p.value
## <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 elevation (Intercept) 78.9 1.43 55.4 0.
## 2 elevation measure -0.104 0.00200 -51.8 0.
## 3 NDVI (Intercept) 4.64 0.388 12.0 2.82e-32
## 4 NDVI measure 1.87 1.70 1.10 2.72e- 1
## 5 shrub_cover (Intercept) 4.94 0.0517 95.5 0.
## 6 shrub_cover measure 5.38 1.54 3.50 4.68e- 4
## 7 slope (Intercept) 5.27 0.0760 69.3 0.
## 8 slope measure -0.0672 0.0223 -3.02 2.54e- 3
## 9 solar (Intercept) 165. 10.9 15.2 1.84e-50
## 10 solar measure -0.000422 0.0000286 -14.7 1.29e-47